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Abstract 

Sex chronnosome dosage connpensation balances honnogametic sex chromosonne expression with autosonnal expression in the 
heterogannetic sex, leading to sex chromosome expression parity between the sexes. If compensation is incomplete, this can lead 
to expression imbalance and sex-biased gene expression. Recent work has uncovered an intriguing and variable pattern of dosage 
compensation across species that includes a lack of complete dosage compensation in ZW species compared with XY species. This has 
led to the hypothesis that ZW species do not require complete compensation or that complete compensation would negatively affect 
their fitness. To date, only one study, a study of the moth Bombyxmori, has discovered evidence for complete dosage compensation in 
a ZW species. We examined another moth species, Manduca sexta, using high-throughput sequencing to survey gene expression in 
the head tissue of males and females. We found dosage compensation to be complete in M. sexta with average expression between 
the Z chromosome in males and females being equal . When genes expressed at very low levels are removed by filtering, we found that 
average autosome expression was highly similar to average Z expression, suggesting that the majority of genes in M. sexta are 
completely dosage compensated. Further, this compensation was accompanied by sex-specific gene expression associated with 
important sexually dimorphic traits. We suggest that complete dosage compensation in ZW species might be more common than 
previously appreciated and linked to additional selective processes, such as sexual selection. More ZW and lepidopteran species should 
now be examined in a phylogenetic framework, to understand the evolution of dosage compensation. 

Key words: Lepidoptera, dosage compensation, sex-biased gene expression, phototransduction, olfaction, mushroom body. 



Introduction 

Sex chromosome dosage compensation is a genetic regulatory 
mechanism that equalizes the expression of the homogametic 
sex chromosome and autosomal genes in the heterogametic 
sex of diploid organisms with a chromosomal-based sex de- 
termination system (Ohno 1967; Mank et al. 2011). In the 
heterogametic sex (males in the XY system and females in 
the ZW system), sex-linked genes on the homogametic sex 
chromosome (X or Z) are present as single copies compared 
with two copies in the homogametic sex (XX or ZZ). In many 
animal species, this has resulted in the evolution of complete 
compensation, a global mechanism of expression regulation 
applied to entire sex chromosomes to balance the expression 
of X/Z and autosomal genes, and which leads to the equali- 
zation of homogametic sex chromosome expression between 
the sexes. For example, in some mammals (including humans). 



gene dose is equalized through the inactivation of one X chro- 
mosome in homogametic females (XX) during early develop- 
ment, coupled with an up-regulation of X-linked genes, or a 
down-regulation of autosomal genes that are functionally 
linked to X-linked genes (Lyon 1961, 1999; Nguyen and 
Disteche 2006; Lin et al. 2007; Johnston et al. 2008; Julien 
et al. 2012). This equalizing of expression level leads to bal- 
anced transcription rates between males and females and 
preserves the integrity of gene networks that include both 
autosomal and sex-linked genes. However, not all mecha- 
nisms of dosage compensation are sex chromosome wide, 
and some species appear to compensate only a subset of 
genes on the X or Z (incomplete dosage compensation; 
Mank 2013). Emerging evidence suggests that patterns of 
dosage compensation are highly variable across sex-determi- 
nation system and species (Mank et al. 201 1). 
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The evolution of sex chronriosonnes proceeds with the dif- 
ferentiation of initially homologous chromosomes that accu- 
mulate antagonistic mutations, leading to selection for 
reduced recombination in the heterogametic sex and the 
gradual degradation of one homolog (Y and W in the XY 
and ZW systems, respectively; Bull 1983; Rice 1987; 
Charlesworth 1991; Charlesworth B and Charlesworth D 
2000). Dosage compensation is thought to have evolved to 
compensate for this degradation (Ohno 1967). There is in- 
creasing evidence for a spectrum in the extent of sex chromo- 
some dosage compensation across species ranging from 
complete to incomplete compensation, applied through 
global and gene-by-gene mechanisms of gene expression reg- 
ulation (Vicoso and Bachtrog 2009). For example, although 
humans display complete dosage compensation, recent evi- 
dence suggests that approximately 15% of genes escape 
compensation on the X chromosome (Carrel and Willard 
2005). Egg laying monotreme mammals (platypus and 
echidna, XY species) show a gene-by-gene pattern of incom- 
plete dosage compensation of the X chromosome, whereas 
marsupials demonstrate X-inactivation and up-regulation of X 
expression through gene-by-gene effects that vary between 
tissues and species (McKay et al. 1987; Cooper et al. 1993; 
Deakin et al. 2009; Julien et al. 2012). Dosage compensation 
in Drosophila (XY) species is complete and achieved through 
the doubling of X chromosome expression in males 
(Mukherjee and Beermann 1965; Belote and Lucchesi 1980; 
Straub et al. 2005). In contrast, the ZW trematode parasite, 
Schistosoma mansoni, has no global mechanism of chromo- 
somal dosage compensation (Vicoso and Bachtrog 2011), 
demonstrating gene-by-gene effects. ZW chickens and 
zebra finches also lack a global mechanism of dosage com- 
pensation (Mank and Ellegren 2009; Zhang et al. 201 0; Ayers 
et al. 2013), yet their patterns of gene-by-gene effects differ 
(Itoh et al. 2010). Further, the flour beetle, Tribolium casta- 
neum (XY), is thought to exhibit global X compensation to 
balance X and autosomal expression, yet female X expression 
is also upregulated leading to sex-biased gene expression 
(Prince et al. 2010). The variability in the degree of dosage 
compensation across species with differentiated sex chromo- 
somes suggests that dosage compensation may not be a re- 
quirement for sex chromosome evolution and that other 
selective forces may shape dosage compensation evolution 
(Mank et al. 201 1 ; Lin et al. 201 2). Additionally, the complex 
phylogenetic distribution of dosage compensation indicates 
that it has evolved independently a number of times and high- 
lights the importance of gene dose regulation to organismal 
fitness (Straub and Becker 2007; Julien et al. 2012; Livernois 
etal. 2013). 

Incomplete dosage compensation produces an expression 
imbalance between males and females and is likely to be re- 
sponsible for a greater proportion of sex-biased expression of 
Z/X-linked genes than previously assumed (Parsch and Ellegren 
2013). Dosage compensation evolved during the evolution of 



sex chromosomes. Thus, the evolution of sex chromosomes, 
dosage compensation, and sex-biased gene expression are 
likely to be constrained by the selective processes acting on 
each. Several factors might influence the evolution of dosage 
compensation and sex-biased gene expression, such as the 
concentration of genes on a chromosome that require 
dosage compensation and the strength of divergent selection 
between the sexes (Livernois et al. 2012). It is unclear how 
often sex-biased gene expression is due to incomplete dosage 
compensation, and although sex is determined by the sex 
chromosomes, sex-specific traits are not necessarily controlled 
by sex chromosome genes (Mank 2009). However, the enrich- 
ment of sex-biased genes on the homogametic sex chromo- 
some is seen frequently in ZW species where dosage 
compensation is rarely complete (Kaiser and Ellegren 2006; 
Arunkumar et al. 2009; Zhang et al. 2010). 

The majority of ZW species demonstrate a pattern of in- 
complete dosage compensation, suggesting that complete 
compensation might be limited to male heterogametic species 
(Mank 2013). Numerous studies, focusing mainly on birds but 
also including snakes and Lepidoptera, have found incomplete 
dosage compensation across ZW species (Mank and Ellegren 
2009; Itoh et al. 2010; Naurin et al. 2012; Adolfsson and 
Ellegren 2013; Vicoso, Emerson, et al. 2013), leading to the 
conclusion that such mechanisms may not exist in the ZW 
system (Wolf and Bryk 2011). This suggests that heteroga- 
metic females might be less sensitive to gene dose or that 
the presence of male-adapted genes on the Z, due to higher 
Z mutation rates and effective population sizes in males (de- 
pendent on several demographic and biological factors, see 
Mank et al. 2010), requires mitigation of antagonistic effects 
through lower female expression (Naurin et al. 2010). This 
means that incomplete dosage compensation might be a 
common cause of sex-biased gene expression in the ZW 
system. It should be noted that incomplete compensation 
leading to homogametic-biased expression is not unique to 
the ZW system. For example, the stickleback, Gasterosteus 
aculeatus (XY), has no chromosome-wide method of compen- 
sation and a higher number of female-biased genes on the X 
chromosome (Leder et al. 2010). Lepidoptera are female het- 
erogametic, and the Indian meal moth, Plodia interpunctella, 
lacks complete Z chromosome dosage compensation 
(Harrison et al. 2012). Incomplete dosage compensation was 
also discovered in the silkworm Bombyx mori using microarray 
data (Zha et al. 2009). However, that result was recently chal- 
lenged when complete Z dosage compensation was shown 
between males and females using high-throughput mRNA 
sequencing (Walters and Hardcastle 2011). The latter study 
is the only example of complete and globally applied dosage 
compensation of a Z chromosome. Walters and Hardcastle 
(201 1 ) also saw that autosomal expression was elevated com- 
pared with Z-linked expression, an unexpected result under 
the expectation of complete compensation. Previous microar- 
ray-based studies of dosage compensation have been called 
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into question (e.g., Xiong et al. 2010), and recent work has 
examined how different methodological aspects of evaluating 
gene expression can affect the conclusions drawn from 
dosage compensation studies (Jue et al. 2013). It is becoming 
clear that best practices need to be developed to examine the 
evolution of dosage compensation in a comparative manner 
(Mank 2013). 

Here, we assess dosage compensation and sex-biased gene 
expression using high-throughput sequencing of mRNA from 
the heads of adult males and females of Manduca sexta, a 
moth species without a fully sequenced genome. We find that 
average gene expression levels between males and females, 
and between the Z chromosome and autosomes, were equal, 
providing evidence for a global mechanism of complete 
dosage compensation along the M. sexta Z chromosome. 
This result highlights the substantial variability in presence 
and degree of dosage compensation in the ZW system in 
general, and Lepidoptera in particular. Further, we found 
that 1,385 contigs were significantly differentially expressed 
(DE) between males and females, with two thirds of these 
genes being female biased. These sex-biased genes were 
distributed throughout the M. sexta genome, with no over- 
representation on the Z chromosome, confirming a lack of 
male-biased Z chromosome expression in M. sexta that is nor- 
mally associated with incomplete dosage compensation in ZW 
species. A functional analysis of sex-biased genes demon- 
strated significant enrichment for a number of sex-related 
functions, including alternatively spliced genes, adult behavior, 
developmental processes, and, interestingly, sensory percep- 
tion of light and neurological development relating to 
olfaction. 

Materials and Methods 

Sequencing and Assembly of the M. sexta Transcriptome 

Manduca sexta larvae were raised on a wheat germ-based 
artificial diet in a climate-controlled chamber at 26 °C, light 
1 6 h/25 °C dark 8 h cycle. To synchronize the eclosion time of 
the adults, the pupae were kept in a climate-controlled cham- 
ber with the following cycles: 14°C, light 16h/14°C dark 8 h 
cycle for 3 weeks, and then at 27 °C, light 12 h/25 °C dark 
1 2 h cycle for 2.5-3 weeks. All pupae were reared separately, 
based on gender, and four replicate samples were collected 
for each sex. Adult heads were collected between 1 5 and 20 h 
after eclosion, then immediately stored in RNAIater at -20 °C 
for RNA extraction. RNA was isolated using an RNeasy Mini Kit 
(Qiagen), following the manufacturer's protocol. The RNA-seq 
libraries were constructed using a strand-specific RNA-Seq 
method according to Zhong et al. (201 1). Briefly, polyA-RNA 
was isolated using Dynabeads Oligo(dT)25 (Invitrogen), then 
simultaneously eluted and fragmented in Superscript III buffer. 
First-strand cDNA was synthesized using Superscript III 
(Invitrogen). Second-strand cDNA was synthesized using 



RNase H (NEB) and DNA polymerase I (NEB) with a dUTP 
mix. The double-stranded cDNA fragments were then sequen- 
tially subjected to end repair, dA tailing, Y-shape adapter liga- 
tion, and Uracil DNA Glycosylase (NEB) treatment. Finally, the 
product was purified and polymerase chain reaction (PGR) 
amplified with indexed PGR primers. Sequencing (50 bp 
single-end reads) was performed on the lllumina HiSeq2000 
platform at Weill Gornell Medical Gollege. 

The eight RNA-seq samples were used to perform a de novo 
transcriptome assembly using the Trinity pipeline (Grabherr 
et al. 2011; Haas et al. 2013). Trinity assembles sequenced 
reads based on sequence similarity using de Brujin graphs, 
constructing reads into contigs that represent genes and alter- 
natively spliced transcripts. Once assembled, each library was 
separately mapped back to the assembly using RSEM (Li and 
Dewey 201 1), and count levels of reads within each genomic 
feature (contigs) were extracted, lllumina reads for each library 
are available as fastq files in the ArrayExpress database (www. 
ebi.ac.uk/arrayexpress, last accessed March 6, 2014) under ac- 
cession number E-MTAB-2066. The FASTA file containing con- 
tigs from the Trinity assembly was deposited in Dryad under 
data identifier doi:10.5061/dryad.gb135. 

Physical Locations of Genes and Dosage Gompensation 

To determine the physical positions of M. sexta genes, de novo 
assembled contigs were compared with the annotated 
genome of the silkworm, B. mori, using the Basic Local 
Alignment Search Tool (Blast; Altschul et al. 1990). There is 
a high level of synteny between the Z chromosomes of but- 
terfly and moth species thus far sequenced {Bombyx, Mita 
et al. 2004; Xia et al. 2004; Danaus, Zhan et al. 2011; 
Heliconius, Dasmahapatra et al. 2012), making it possible to 
assign the putative physical locations of M. sexta Z and auto- 
somal contigs through orthology to B. mori (Harrison et al. 
2012). A reciprocal best-hit approach was used to determine 
highly similar 1:1 orthologs using a strict e-value cut off of 
1 X 10~^°. Manduca sexta contigs were assigned putative 
chromosome locations from their corresponding 1:1 orthologs 
in B. moh. All contigs analyzed were more than 200 bp in 
length, and fragments per kilobase per million mapped 
reads (FPKM), a conversion factor that accounts for contig 
length when considering expression levels, was calculated 
for each contig. Libraries were normalized using the trimmed 
mean of A// values (TMM) normalization method in the NOISeq 
R package (Tarazona et al. 2011) before assessing dosage 
compensation. The four replicate male and female libraries 
were tested for uniformity using a Spearman's rank order cor- 
relation test. FPKM values for same-sex replicate libraries 
were highly correlated (pairwise Spearman's p > 0.93; 
P<2.2x10~^^ between female replicates, and p>0.91; 
P<2.2 X 10~^^ between males), thus male and female repli- 
cates could be reliably combined to compare expression levels 
across Z and autosomal contigs. From the Blast results, contigs 
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were partitioned into test sets for putative M. sexta autosonnal 
and Z chromosome genes. Male:female gene expression ratio 
was investigated by log2 transforming FPKM for all contigs in 
the test set, as well as for genes with 1:1 orthologs in B. mori, 
to ensure there was no bias caused by removing contigs with 
reduced orthology. 

Analyses of chromosome-wide expression data have been 
criticized for the potential bias introduced by a lack of robust 
filtering (Xiong et al. 2010; Kharchenko et al. 201 1). Filtering 
of microarray data is common due to levels of background 
noise, below which gene expression cannot be accurately 
called. RNA sequencing can also be biased, depending on 
factors such as sequencing depth, and can influence conclu- 
sions as to the extent of dosage compensation (Kharchenko 
et al. 2011; Walters and Hardcastle 2011; Jue et al. 2013). 
Therefore, we explored our data using two different filtering 
methods used in previous investigations of dosage compen- 
sation, as well as assessing the unfiltered data set. First, filter- 
ing was performed as per Jue et al. (2013) to remove outlier 
contigs from the distribution of expression levels. The FPKM of 
each contig was log2 transformed, and outliers were removed 
if they were either 1) more than 1 .5 times the mid-50 percen- 
tile above the 75^^ percentile or 2) the same distance below 
the 25^^ percentile, thus producing data that approximates a 
normal distribution. Second, our test set of contigs was fil- 
tered (using unlogged FPKM values) to remove low-expressed 
genes as per Harrison et al. (2012). We explored three sepa- 
rate data sets with different stringencies of this filtering 
method, discarding contigs that had expression levels below 
2, 3, and 4 FPKM in at least four groups. Using these differing 
filtering methods for contigs with putative physical locations, 
we compared chromosome-wide expression patterns of the Z 
and autosomes, between males and females. Significance 
testing of distributions of FPKM values was performed using 
a two-sample Wilcoxon rank sum test, and for normally dis- 
tributed log2 transformed data, Mests were employed. 

Gene Expression Analysis and Functional Enrichment 

Raw count values from read mapping to the de novo assembly 
were extracted using Trinity perl scripts for each of the four 
male and four female libraries. The edgeR package in R was 
used to perform differential expression analysis to identify DE 
contigs between males and females (Robinson et al. 2010). 
Normalization was performed in edgeR, using the TMM 
method and included adjustments for library size to account 
for between-sample variation in sequencing depth. For data 
sets with replicate samples, edgeR performs a generalized 
linear model using a negative binomial distribution to identify 
DE genes. Significantly, DE contigs were determined by using 
the false discovery rate (FDR) of Storey and Tibshirani (2003) to 
adjust gene-specific P values for multiple tests, taking a signif- 
icance threshold of 5% FDR. A heatmap was constructed 
using the normalized counts of significant contigs to plot 



changes in moderated log2 counts per million (log CPM) 
across libraries, with the R package Heatplus (Ploner 2012). 
To address the possibility that edgeR does not sufficiently nor- 
malize libraries that have unequal sequencing depth, we per- 
formed differential analysis on a subsampled data set. Briefly, 
we down-sampled sequence libraries to the smallest library 
size, thereby equalizing all library sizes. Down-sampled librar- 
ies were mapped to the assembled contigs and differential 
expression analysis was repeated, as for the full data set. 

To examine the role of coexpressed genes, functional en- 
richment tests can provide a statistically based summary of 
functional terms that occur more often than chance in a 
gene test set. Functional enrichment was performed using 
orthologous genes in Drosophila melanogaster, obtained 
through a Blast comparison of all M. sexta contigs to the 
D. melanogaster genome (assembly BDGP5, release 71). A 
Blast hit was considered orthologous when the best top hit 
passed an e-value threshold of less than 1 x 10~^. Functional 
enrichment of these orthologs was performed in GOrilla (Eden 
et al. 2009) and DAVID v6.7 (Dennis et al. 2003; Huang et al. 
2009), testing the functions of Drosophila orthologs of signif- 
icant M. sexta sex-biased contigs against the D. melanogaster 
gene set. The GOrilla functional enrichment test produces a 
list of significant individual gene ontology (GO) terms, 
whereas DAVID performs enrichment tests using a number 
of different functional terms, including GO terms, SwissProt 
keywords, and UniProt sequence features. Additionally, 
DAVID performs a fuzzy heuristic partitioning that groups 
GO terms by their degree of overlapping gene sets, reducing 
the presence of repetitive or nested terms. GOrilla GO terms 
and DAVID clusters were considered significant with an 
FDR<5% and with an enrichment score (the geometric 
mean of annotation P values) of more than 1 .3, respectively. 
DAVID clusters were filtered by hand to remove obvious over- 
lapping clusters. 

Results 

Dosage Compensation 

The M. sexta de novo transcriptome assembly resulted in a 
total of 31,459 gene contigs, with a total combined length of 
23.9 Mb and an average per nucleotide transcriptome cover- 
age of 34.7 X for uniquely mapped reads across eight libraries, 
which ranged in size from 1 6.6 to 24.3 million reads. To assess 
the presence or absence of dosage compensation in M. sexta, 
assembled contigs were Blast searched against the B. mori 
genome to obtain the physical locations of 1:1 orthologs. 
We found 575 putative Z chromosome and 1 3,799 autosomal 
M. sexta contigs that had 1:1 orthology to B. mori genes and 
were more than 200 bp length. Filtering was performed to 
reduce bias from extreme expression levels, using two differ- 
ent methods. The number of retained putative M. sexta Z 
contigs under these different methods ranged from 451 to 
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568, and 11,956 to 13,689 for autosomal contigs (table 1), 
depending on the stringency of filtering. All pairwise correla- 
tions of FPKM values between replicate samples were highly 
significant (P<2.2 x 10~^^), and thus, averaged expression 
levels could be assessed for dosage compensation. 

In general, comparisons of global gene expression between 
males and females revealed minimal differences, and this was 
true for both the Z chromosome and for autosomes. Thus, 
these observations indicate complete dosage compensation in 
M. sexta (figs. 1 and 2). The distributions of male:female gene 
expression ratios between contigs on the Z chromosome and 
the autosomes overlapped completely (fig. la). Considering all 
assembled contigs together (i.e., regardless of physical loca- 
tion) resulted in a unimodal distribution (fig. ^b), indicating 
male:female parity. Autosomal expression differences be- 
tween males and females were small (supplementary table 

51, Supplementary Material online; fig. 2), matching the ex- 
pectation of autosomal parity between males and females. 
However, the average expression level of autosome compared 
with Z chromosome contigs did show a small but vari- 
able difference that was dependent on filtering stringency 
(table 1). 

Average Z chromosome expression between males and fe- 
males was nearly identical and did not differ significantly at 
any level of filtering (Wilcoxon test; P> 0.05 for all compari- 
sons; table 1, supplementary table SI, Supplementary 
Material online), indicating that the presence or level of filter- 
ing stringency did not affect the conclusion that the Z chro- 
mosome shows complete dosage compensation in M. sexta. 
However, a significant, though small, difference was seen be- 
tween average autosome and Z expression in both males and 
females when data were unfiltered or filtered at low strin- 
gency (unfiltered, log2 transformed outlier removal and data 
filtered at 2 FPKM; fig. 2a; table 1 ; supplementary fig. SI and 

52, Supplementary Material online). This pattern of increased 
autosomal expression was strongest in the unfiltered data set 



(Wilcoxon tests; female P= 1 .52 x 1 0~^; male 
P= 5.78 X 10"^), but average expression differences and 
strength of significance reduced rapidly with a small increase 
in filtering threshold: Although differences were significant 
when filtering at 2 FPKM (female P= 6.26x10"^; male 
P=2.02 X 10~^), they were only marginally significant at 3 
FPKM (female P=0.01; male P=0.05) and nonsignificant 
when filtering at 4 FPKM (female P=0.12; male P=0.96; 
fig. 2b: table 1; supplementary figs. S1-S3, Supplementary 
Material online). Filtering at the 4 FPKM threshold resulted 
in removal of 21.6% of Z chromosome contigs and 13.4% 
of autosomal contigs (supplementary table S2, Supplementary 
Material online) and led to Z:autosome parity. Thus, our anal- 
ysis reveals that when data from genes expressed at very low 
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Log2 maleifemale ratio of mean expression level (FPKM) 

Fig. 1. — Histogram showing log2 maleifemale ratio of expression level 
(FPKM) for Manduca sexta contigs located on the Z chromosome (red) and 
the autosomes (blue; a), and all contigs regardless of physical location (b). 



Table 1 

Effect of Filtering Method and Stringency on the Number of Genes in Each Analysis and the Presence of Dosage Compensation 



Filtering Method 


Number of Z 


Number of 


Z, M:F 


M, Z:A 


F, Z:A Parity*** 




Genes Retained 


Autosomal Genes 


Parity* 


Parity** 








Retained 








Unfiltered 


575 


13,799 


0.81 


5.78 X 10"°^ 


1.52 X 10"°^ 


outlier removal 


568 


13,689 


0.80 


2.26 X 10-°"^ 


3.06 X 10"°^ 


Remove genes < 2 FPKM 


559 


13,551 


0.83 


2.02 X 10"°"^ 


6.26 X 10"°^ 


Remove genes < 3 FPKM 


503 


12,874 


0.93 


0.05 


0.01 


Remove genes < 4 FPKM 


451 


11,956 


0.96 


0.96 


0.12 



Note. — Filtering methods are outlined in the Materials and Methods section, genes were assigned location using 1:1 orthologs in Bombyx mori. Sex chromosome dosage 
compensation was assessed by comparing Z chromosome expression between males and females, and Z expression was compared with autosome expression for males and 
females separately. All significance tests were Wilcoxon tests except for the normally distributed outlier method of filtering, where a t-test was employed. Significance 
indicates disparity in expression levels. 

*P values for Z, M:F parity compare average Z chromosome expression between males and females. 

**P values for M, Z:A parity compare average male expression between Z and autosomal genes. 

***P values for F, Z:A parity compare average female expression between Z and autosomal genes. 
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Fig. 2. — ^Average Z-linked and autosomal contig expression levels (FPKM) across four replicate males and females using the log2 outlier removal method 
of filtering (a) and removal of contigs with less than 4 FPKM (b). Black lines are the median of the FPKM distribution across contigs, boxes show the 
interquartile range, whiskers extend to 1 .5 x the interquartile range and notches approximate the 95% confidence intervals of the medians. Overlapping 
notches are evidence for the similarity of median values. 



levels are rennoved by filtering, the remaining majority are fully 
dosage compensated and equal between the sexes. The out- 
lier removal filtering method produced a test set of 568 Z 
chromosome contigs and 13,689 autosome contigs for com- 
parison, which had a mean log2 FPKM of 3.65, 3.35, 3.6, and 
3.33 for autosome female, Z female, autosome male, and Z 
male, respectively (fig. 2a, supplementary table S1, 
Supplementary Material online). Median expression levels 
when filtering at a threshold of 4 FPKM were 12.75, 12.55, 
12.7, and 12.52 for autosome female, Z female, autosome 
male, and Z male, respectively (fig. 2b; supplementary table 
SI, Supplementary Material online). 

Differential Expression and Functional Analyses 

Through differential expression analysis of males and females, 
we identified a total of 1,385 sex-biased contigs from the full 
data set (FDR < 5%). The majority of these sex-biased contigs 
(960) were upregulated in females (female-biased; fig. 3a). 
Because unequal library sizes may influence differential ex- 
pression results using edgeR, we also performed differential 
expression on down-sampled libraries, which were of equal 
size, to compare to the full data set results. Analysis of this 
subsampled data demonstrated very similar results to that of 
the full data set. In total, 1,310 genes were significantly DE 
(FDR < 0.05), with 906 being female biased and 404 male 
biased, overlapping with the full data set results considerably 
(supplementary fig. S4, Supplementary Material online). 
Mapping efficiencies of the subsampled data were high and 
equal between libraries (>91 % for all libraries; supplementary 
table S3, Supplementary Material online). Therefore, ineffi- 
cient mapping and unequal library sizes are unlikely to have 



contributed to our full data set differential expression results. 
Thus, the full data set results were used for all further analyses. 

To examine the distribution of sex-biased genes across the 
genome, DE contigs from the full data set (1,385 contigs) 
were assigned a putative physical location when they had 
1:1 orthology to B. mori. Of the 1,385 sex-biased contigs, 
604 had 1:1 orthologs with positional information. The phys- 
ical locations of sex-biased contigs did not exhibit a bias to- 
ward the Z chromosome, as would be expected in the 
absence of dosage compensation, and sex-biased genes 
were distributed throughout the genome (fig. 3b). Average 
expression levels of sex-biased contigs between males and 
females from the Z chromosome and autosomes were also 
examined (570 autosomal contigs and 34 Z-linked contigs; 
supplementary fig. S5, Supplementary Material online). The 
average expression of sex-biased contigs from the Z chromo- 
some appeared higher than that of contigs from the auto- 
somes; however, this effect was not significant (Wilcoxon 
test; P>0.05). 

Functional enrichment of sex-biased contigs was per- 
formed using D. melanogaster orthologs obtained from a 
Blast search of D. melanogaster cDMAs. In total, 5,793 contigs 
had good hits to genes in the D. melanogaster genome, rep- 
resenting 18.4% of M. sexta contigs. Of these 5,793 contigs, 
1 91 were sex biased in the differential expression analysis, and 
thus, functional information from their D. melanogaster 
orthologs could be used to perform functional enrichment 
tests. Generally, both DAVID and GOrilla demonstrated signif- 
icant enrichment of similar GO terms, with DAVID producing a 
more succinct set of functional clusters (table 2). 

Nineteen DAVID clusters demonstrated an enrichment 
score of more than 1.3, 14 of which had nonredundant 
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Fig. 3. — Heatmap of significantly DE genes (FDR < 0.05) plotted using log2 counts per million (log CPM) per gene for each Manduca sexta library (males 
1-4 and females 1-4, as labeled; a). Color bar legend indicates the degree of log CPM change between males and females, (b) Physical locations of sex- 
biased genes determined using 1:1 orthologs in Bombyx mori. Chromosome numbers are shown, with the Z chromosome being chromosome 1, as 
indicated. The size of each chromosome segment in the pie chart represents the relative abundance of sex-biased genes, normalized to chromosome size. 



functions (table 2). Many of these GO terms were develop- 
mental (e.g., GO:0060284, regulation of cell development); 
however, several involved the sensing and parsing of external 
cues and changes in behavior. These included GO terms for 
the response to an external or abiotic stimulus that were spe- 
cifically involved in sensory perception of light stimulus 
(GO:0050953; e.g., ninaC, inaD, and trpi). Behavior-related 
GO terms included behavior itself (GO:0007610; e.g., cpo, 
b, and Syn) and adult locomotory behavior (GO:0008344; 
e.g., tuti, cac, and t). Other clusters included annotations in- 
volved in the regulation of development, specifically regula- 
tion of nervous system development (GO:0051960; e.g., 
CadN, shot, and sdk), and neurological development that in- 
cluded the term mushroom body development (GO:0016319; 
e.g., RhoGAPp190, start, and chinmo). The regulation of sex- 
specific gene expression was represented by two clusters, one 
containing annotations for alternative splicing and the other 
for RNA editing. GOrilla produced 1 02 significant GO terms in 
total (FDR<5%; supplementary table S4, Supplementary 
Material online). 

To visualize gene expression changes for contigs of a spe- 
cific function, we examined the expression levels of four func- 
tional groups of contigs, chosen based on their inclusion in 
significant GO terms in the DAVID functional enrichment test. 
The four GO terms were phototransduction (8 genes; 
GO:0007602), deactivation of rhodopsin mediated signaling 
(5 genes; GO:0016059), regulation of cell development (11 
genes; GO:0060284), and behavior (17 genes; GO:0007610). 
Mean expression levels (FPKM) of genes within each GO term 
were calculated across four replicates. Genes within each 



functional annotation were all upregulated in females com- 
pared with males, supporting the general pattern of female- 
biased gene expression (supplementary fig. S6 and table S5, 
Supplementary Material online). 

Discussion 

Complete sex chromosome dosage compensation has not 
been unequivocally observed in a species with a ZW sex de- 
termination system, leading to the supposition that complete 
compensation of sex chromosome expression only occurs in 
male heterogametic species. In fact, only one study has shown 
complete dosage compensation in a ZW species, following a 
reassessment of previously published data (Zha et al. 2009; 
Walters and Hardcastle 201 1 ). Our results show that complete 
dosage compensation has evolved in the ZW species M. sexta, 
alongside sex-specific gene expression. Regardless of the 
method and stringency of data filtering performed, the aver- 
age expression of Z-linked contigs between males and females 
was essentially identical, suggesting the complete compensa- 
tion of Z-linked genes. Further, in agreement with Walters and 
Hardcastle (201 1), we found the average expression level of 
autosomal contigs to be higher than that of Z-linked contigs. 
However, this imbalance quickly disappears when contigs 
with low expression levels are removed (supplementary fig. 
SI -S3, Supplementary Material online, vs. fig. 2b), suggesting 
that the majority of genes in the M. sexta genome are fully 
dosage compensated between the Z chromosome and 
autosomes. 
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DAVID Functional Enrichment Results 
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Note. — All GO term clusters presented have an enrichment score of more than 1.3 (enrichment scores are shown underneath the cluster number). Enrichment scores are 
the geometric mean of GO term P values within a cluster. Each cluster is presented as the top three representative GO terms for that cluster and includes the number of sex- 
biased genes, annotated with each functional GO term, the P value and FDR. Annotation term notation: GOTERM, GO term; BP, biological process; CC, cellular component; 
MF, molecular function; FAT; DAVID database filtered for specific GO terms; SP_PIR_KEYWORDS, Swiss-Prot Protein Information Resource Keywords; UP_SEQ_FEATURE, 
UniProt sequence feature. 
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Dosage Compensation 

Dosage compensation has previously been considered as nec- 
essary for sex chromosome evolution, to equalize changes in 
gene product concentration as one chromosome degener- 
ates. However, dosage compensation studies are often con- 
tradictory (Livernois et al. 2013), and patterns of dosage 
compensation are varied across species (fig. 4), ranging from 
incomplete compensation to the hyperexpression of sex chro- 
mosomes in both sexes (Prince et al. 2010). Incomplete 
dosage compensation exists in XY species, but is by far 
more common in the ZW system. It has been suggested 
that ZW species have not evolved complete compensation 
because either 1) heterogametic females are not highly ad- 
versely affected by unequal dose and require fewer genes on 
the Z to be compensated or 2) the presence of male-adapted 
genes on the Z has led to a trade-off between maintaining 
gene network integrity and expressing maladaptive genes in 
females (Mank and Ellegren 2009; Naurin et al. 2010; Mank 
et al. 201 1 ; Naurin et al. 201 2). However, we have discovered 
complete dosage compensation in a ZW species, alongside 
genome-wide sex-biased gene expression. This suggests that 
complete dosage compensation is required to balance the 
expression of key gene networks throughout the genome of 
M. sexta. The evolution of complete compensation may have 
allowed for the evolution of sex-biased expression through 
gene-by-gene mechanisms driven by sex-specific selective 
pressures such as sexual selection or sex-specific natural selec- 
tion. Alternatively, sex-biased expression might have evolved 
alongside dosage compensation, under sexually antagonistic 
selection (Vicoso, Kaiser, et al. 2013). This might occur when 



specific genome-wide loci are under strong selection for op- 
timal expression levels between the sexes, with dosage com- 
pensation of Z-linked genes allowing the fine tuning of sex- 
specific dose. The presence of complete compensation in a 
related species, B. mori, indicates that complete dosage com- 
pensation might have evolved in a common ancestor of these 
two species. Thus, further studies are necessary to determine 
whether complete compensation is common in moths, and in 
ZW species generally, and to improve our understanding of 
the phylogenetic distribution of compensation and the mech- 
anisms that underlie it (fig. 4). 

The resolution of conflict during sex chromosome evolution 
through gene-by-gene regulation of sex-linked genes is 
thought to be a counter mechanism to sex chromosome dif- 
ferentiation, allowing for sex determination with homomor- 
phic sex chromosomes (e.g., in the Emu; Vicoso, Kaiser, et al. 
2013) and thus not requiring the evolution of dosage com- 
pensation (e.g.. Boa constrictor; Vicoso, Emerson, et al. 201 3). 
This is an interesting aspect of the dynamics of sex chromo- 
some evolution wherein both dosage compensation and sex- 
biased gene expression are essentially the result of sexual con- 
flict over optimal trait values for males and females, yet the 
target loci of these processes may differ (Mank et al. 201 1). 
Therefore, complete sex chromosome dosage compensation 
may not be required during sex chromosome evolution in ZW 
species but in some instances may have evolved alongside 
gene-by-gene mechanisms for key sex-specific loci across 
the genome. The signal of female-biased gene expression 
seen here is unusual for ZW species, which generally demon- 
strate expression of male-biased Z-linked genes due to 
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Fig. 4. — Putative models of X/Z chromosomal dosage compensation for representative XY and ZW species. Chromosomes represent the X/Z (large 
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inconnplete dosage connpensation of the Z chronnosonne (e.g., 
in birds; Naurin et al. 2012). Thus, female-biased gene expres- 
sion is nnore connnnon in XY species where dosage connpen- 
sation is more frequently complete. It is likely that sex-biased 
gene expression is a result of sexual selection in M. sexta, as 
sex-biased genes were found to function in known sexually 
dimorphic traits (see later). This indicates that some patterns of 
sex-specific expression seen on the Z in other ZW species may 
also be due to sex-specific selection (Mank et al. 201 1). 

Sex-Biased Genes in M. sexta 

We found a number of contigs to be DE between the sexes, 
and enrichment analyses suggested that these differences 
were associated with a number of sexually dimorphic func- 
tions. Manduca sexta is sexually dimorphic in antennal lobe 
morphology, with sex-specific differences in size and position 
of glomeruli (olfactory centers of the insect antennal lobe; 
Rospars and Hildebrand 2000). This underlies a general dimor- 
phism in olfactory system and sensory physiology in M. sexta 
(Shields and Hildebrand 2001). These differences correlate 
with our expression data from male and female heads, with 
significant enrichment of genes relating to the regulation of 
nervous system and mushroom body development. 
Mushroom bodies are the next step in olfaction pathways 
after the antennal lobe, and they translate olfactory sensory 
information into learned behavioral responses (Matsumoto 
and Hildebrand 1981; Daly and Smith 2000; Zars 2000; 
Caron et al. 2013). Sex-biased genes in M. sexta were func- 
tionally enriched for sensory perception and behavior, includ- 
ing five genes involved in the sensory perception of smell 
(GO:0007608; rdgA, rdgB, eag, norpA, and IrSA). 
Interestingly, there was a strong signal for genes involved in 
visual perception, suggesting that vision might also be an im- 
portant sexually dimorphic trait in this species. Manduca sexta 
is a crepuscular/nocturnal species and thus vision might be of 
particular importance to female-specific behavior in low light 
environments, for example, in the detection of suitable host 
plants for oviposition. Indeed, the majority of sex-biased gene 
functions were toward female expression, suggesting that fe- 
males have evolved specific expression of genes involved in 
visual and olfactory sensory perception and adult behavior. 
Gene expression has been previously surveyed in the antennae 
of M. sexta adults using microarray data (Grosse-Wilde et al. 
201 1). Grosse-Wilde et al. (201 1) discovered a female bias in 
the number of transcripts found in the antennal tissue with 
729 transcripts present in the female compared with 348 in 
the male. Together these data suggest that sex-specific selec- 
tion has acted to regulate genes across the genome related to 
sensory perception and behavior in M. sexta. 

Conclusion 

We examined dosage compensation using two methods of 
filtering and several different stringencies for filtering contigs 



with low expression levels. Increased stringency in filtering led 
to increased parity of Z:autosome expression. Walters and 
Hardcastle (2011) also saw disparity between Z-linked and 
autosome expression; however, we have shown that this is 
only true for a minority of genes. Thus, filtering and analysis 
methods are important for a thorough investigation of the 
extent of dosage compensation. Analytical protocols including 
the initial experimental design, for example, the developmen- 
tal stage of the organism and the tissue assayed, and the 
sequencing and mapping or assembly method can all influ- 
ence conclusions on dosage compensation (Jue et al. 2013). 
Our result of complete dosage compensation between male 
and female Z-linked genes held true across all filtering meth- 
ods and was robust to potential errors in determining physical 
locations (i.e., the unimodal distribution of all male:female 
contigs; fig. lb). Further, complete chromosome and sex 
parity was true for at least three-quarters of the most highly 
expressed genes. The sampling of adult heads allowed us to 
survey the effect of dosage compensation after sex determi- 
nation and assay sex-specific gene expression at an important 
life stage, which included the sexually dimorphic antennal 
tissue. Only a subset of orthologous sex-biased contigs was 
assessed for function, due to few Blast hits to Drosophila 
genes. Ascertaining gene functions for nonmodel organisms 
can often be a challenge when comparing genes with that of 
a divergent model species. Thus, although we recovered a 
number of sex-biased genes with functional annotations for 
known sexually dimorphic traits in M. sexta, it is likely that 
some important functions were missed. The improved anno- 
tation of lepidopteran genomes, and the increasing number of 
genomes being sequenced, will aid in the improvement of 
functional annotation of M. sexta genes. 

The evolution of sex chromosomes, dosage compensation, 
sex-specific gene expression, and sex determination all involve 
the specific regulation of sex-related genes. Each process may 
be important at different developmental stages, and the ex- 
pression of these genes as networks requires their tight regu- 
lation across the genome. We find that M. sexta fully 
compensates for Z chromosome dose between males and fe- 
males and with autosomal expression. This indicates that M. 
sexta requires balance in genome-wide gene network expres- 
sion, and this balance might be necessary for the sex-biased 
expression of autosomal genes that underlie important sexu- 
ally dimorphic traits. We propose that sex-specific selection 
may be an important influence on the evolution of complete 
dosage compensation in ZW species and that additional ZW 
taxa should be examined to determine the extent of compen- 
sation and sex-biased expression in a phylogenetic framework. 

Supplementary Material 

Supplementary figures S1-S6 and tables S1-S8 are available 
at Genome Biology and Evolution online (http://www.gbe. 
oxfordjournals.org/). 
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